﻿#region << 版 本 注 释 >>
/*----------------------------------------------------------------
 * 版权所有 (c) 2022 北京超维景生物科技有限公司 保留所有权利。
 * 
 * 创建者：huangyang
 * 电子邮箱：huangyang@tvscope.cn
 * 创建时间：2023/3/13 17:30:27
 * 版本：V1.0.0
 * 描述：
 *
 * ----------------------------------------------------------------
 * 修改人：
 * 时间：
 * 修改说明：
 *
 * 版本：V1.0.1
 *----------------------------------------------------------------*/
#endregion << 版 本 注 释 >>

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace ImageK.Process
{
    /** Autothresholding methods (limited to 256 bin histograms) from the Auto_Threshold plugin 
    (http://fiji.sc/Auto_Threshold) by G.Landini at bham dot ac dot uk). */
    public class AutoThresholder
    {
        private static string[] mStrings;

        public enum Method
        {
            Default,
            Huang,
            Intermodes,
            IsoData,
            IJ_IsoData,
            Li,
            MaxEntropy,
            Mean,
            MinError,
            Minimum,
            Moments,
            Otsu,
            Percentile,
            RenyiEntropy,
            Shanbhag,
            Triangle,
            Yen
        };

        public static string[] getMethods()
        {
            if (mStrings == null)
            {
                var list = new List<string>();
                foreach (int myCode in Enum.GetValues(typeof(Method)))
                {
                    list.Add(Enum.GetName(typeof(Method), myCode));//获取名称
                    // string strVaule = myCode.ToString();//获取值
                }

                mStrings = list.ToArray();
            }
            return mStrings;
        }

        /** Calculates and returns a threshold using the specified
        method and 256 bin histogram. */
        public int getThreshold(Method method, int[] histogram)
        {
            if (histogram == null)
                throw new ArgumentException("Histogram is null");
            if (histogram.Length != 256)
                throw new ArgumentException("Histogram length not 256");
            int threshold = 0;
            switch (method)
            {
                case Method.Default: threshold = defaultIsoData(histogram); break;
                case Method.IJ_IsoData: threshold = IJIsoData(histogram); break;
                case Method.Huang: threshold = Huang(histogram); break;
                case Method.Intermodes: threshold = Intermodes(histogram); break;
                case Method.IsoData: threshold = IsoData(histogram); break;
                case Method.Li: threshold = Li(histogram); break;
                case Method.MaxEntropy: threshold = MaxEntropy(histogram); break;
                case Method.Mean: threshold = Mean(histogram); break;
                case Method.MinError: threshold = MinErrorI(histogram); break;
                case Method.Minimum: threshold = Minimum(histogram); break;
                case Method.Moments: threshold = Moments(histogram); break;
                case Method.Otsu: threshold = Otsu(histogram); break;
                case Method.Percentile: threshold = Percentile(histogram); break;
                case Method.RenyiEntropy: threshold = RenyiEntropy(histogram); break;
                case Method.Shanbhag: threshold = Shanbhag(histogram); break;
                case Method.Triangle: threshold = Triangle(histogram); break;
                case Method.Yen: threshold = Yen(histogram); break;
            }
            if (threshold == -1) threshold = 0;
            return threshold;
        }

        public int getThreshold(string mString, int[] histogram)
        {
            // throws an exception if unknown argument
            int index = mString.IndexOf(" ");
            if (index != -1)
                mString = mString.Substring(0, index);
            Method method = Enum.Parse<Method>(mString); 
            return getThreshold(method, histogram);
        }

		int Huang(int[] data)
		{
			// Implements Huang's fuzzy thresholding method 
			// Uses Shannon's entropy function (one can also use Yager's entropy function) 
			// Huang L.-K. and Wang M.-J.J. (1995) "Image Thresholding by Minimizing  
			// the Measures of Fuzziness" Pattern Recognition, 28(1): 41-51
			// M. Emre Celebi  06.15.2007
			// Ported to ImageJ plugin by G. Landini from E Celebi's fourier_0.8 routines
			int threshold = -1;
			int ih, it;
			int first_bin;
			int last_bin;
			double sum_pix;
			double num_pix;
			double term;
			double ent;  // entropy 
			double min_ent; // min entropy 
			double mu_x;

			/* Determine the first non-zero bin */
			first_bin = 0;
			for (ih = 0; ih < 256; ih++)
			{
				if (data[ih] != 0)
				{
					first_bin = ih;
					break;
				}
			}

			/* Determine the last non-zero bin */
			last_bin = 255;
			for (ih = 255; ih >= first_bin; ih--)
			{
				if (data[ih] != 0)
				{
					last_bin = ih;
					break;
				}
			}
			term = 1.0 / (double)(last_bin - first_bin);
			double[] mu_0 = new double[256];
			sum_pix = num_pix = 0;
			for (ih = first_bin; ih < 256; ih++)
			{
				sum_pix += (double)ih * data[ih];
				num_pix += data[ih];
				/* NUM_PIX cannot be zero ! */
				mu_0[ih] = sum_pix / num_pix;
			}

			double[] mu_1 = new double[256];
			sum_pix = num_pix = 0;
			for (ih = last_bin; ih > 0; ih--)
			{
				sum_pix += (double)ih * data[ih];
				num_pix += data[ih];
				/* NUM_PIX cannot be zero ! */
				mu_1[ih - 1] = sum_pix / (double)num_pix;
			}

			/* Determine the threshold that minimizes the fuzzy entropy */
			threshold = -1;
			min_ent = double.MaxValue;
			for (it = 0; it < 256; it++)
			{
				ent = 0.0;
				for (ih = 0; ih <= it; ih++)
				{
					/* Equation (4) in Ref. 1 */
					mu_x = 1.0 / (1.0 + term * Math.Abs(ih - mu_0[it]));
					if (!((mu_x < 1e-06) || (mu_x > 0.999999)))
					{
						/* Equation (6) & (8) in Ref. 1 */
						ent += data[ih] * (-mu_x * Math.Log(mu_x) - (1.0 - mu_x) * Math.Log(1.0 - mu_x));
					}
				}

				for (ih = it + 1; ih < 256; ih++)
				{
					/* Equation (4) in Ref. 1 */
					mu_x = 1.0 / (1.0 + term * Math.Abs(ih - mu_1[it]));
					if (!((mu_x < 1e-06) || (mu_x > 0.999999)))
					{
						/* Equation (6) & (8) in Ref. 1 */
						ent += data[ih] * (-mu_x * Math.Log(mu_x) - (1.0 - mu_x) * Math.Log(1.0 - mu_x));
					}
				}
				/* No need to divide by NUM_ROWS * NUM_COLS * LOG(2) ! */
				if (ent < min_ent)
				{
					min_ent = ent;
					threshold = it;
				}
			}
			return threshold;
		}

		bool bimodalTest(double[] y)
		{
			int len = y.Length;
			bool b = false;
			int modes = 0;

			for (int k = 1; k < len - 1; k++)
			{
				if (y[k - 1] < y[k] && y[k + 1] < y[k])
				{
					modes++;
					if (modes > 2)
						return false;
				}
			}
			if (modes == 2)
				b = true;
			return b;
		}

		int Intermodes(int[] data)
		{
			// J. M. S. Prewitt and M. L. Mendelsohn, "The analysis of cell images," in
			// Annals of the New York Academy of Sciences, vol. 128, pp. 1035-1053, 1966.
			// ported to ImageJ plugin by G.Landini from Antti Niemisto's Matlab code (GPL)
			// Original Matlab code Copyright (C) 2004 Antti Niemisto
			// See http://www.cs.tut.fi/~ant/histthresh/ for an excellent slide presentation
			// and the original Matlab code.
			//
			// Assumes a bimodal histogram. The histogram needs is smoothed (using a
			// running average of size 3, iteratively) until there are only two local maxima.
			// j and k
			// Threshold t is (j+k)/2.
			// Images with histograms having extremely unequal peaks or a broad and
			// flat valleys are unsuitable for this method.

			int minbin = -1, maxbin = -1;
			for (int i = 0; i < data.Length; i++)
				if (data[i] > 0) maxbin = i;
			for (int i = data.Length - 1; i >= 0; i--)
				if (data[i] > 0) minbin = i;
			int length = (maxbin - minbin) + 1;
			double[] hist = new double[length];
			for (int i = minbin; i <= maxbin; i++)
				hist[i - minbin] = data[i];

			int iter = 0;
			int threshold = -1;
			while (!bimodalTest(hist))
			{
				//smooth with a 3 point running mean filter
				double previous = 0, current = 0, next = hist[0];
				for (int i = 0; i < length - 1; i++)
				{
					previous = current;
					current = next;
					next = hist[i + 1];
					hist[i] = (previous + current + next) / 3;
				}
				hist[length - 1] = (current + next) / 3;
				iter++;
				if (iter > 10000)
				{
					threshold = -1;
					IJ.log("Intermodes Threshold not found after 10000 iterations.");
					return threshold;
				}
			}

			// The threshold is the mean between the two peaks.
			int tt = 0;
			for (int i = 1; i < length - 1; i++)
			{
				if (hist[i - 1] < hist[i] && hist[i + 1] < hist[i])
				{
					tt += i;
					//IJ.log("mode:" +i);
				}
			}
			threshold = (int)Math.Floor(tt / 2.0);
			return threshold + minbin;
		}

		int IsoData(int[] data)
		{
			// Also called intermeans
			// Iterative procedure based on the isodata algorithm [T.W. Ridler, S. Calvard, Picture 
			// thresholding using an iterative selection method, IEEE Trans. System, Man and 
			// Cybernetics, SMC-8 (1978) 630-632.] 
			// The procedure divides the image into objects and background by taking an initial threshold,
			// then the averages of the pixels at or below the threshold and pixels above are computed. 
			// The averages of those two values are computed, the threshold is incremented and the 
			// process is repeated until the threshold is larger than the composite average. That is,
			//  threshold = (average background + average objects)/2
			// The code in ImageJ that implements this function is the getAutoThreshold() method in the ImageProcessor class. 
			//
			// From: Tim Morris (dtm@ap.co.umist.ac.uk)
			// Subject: Re: Thresholding method?
			// posted to sci.image.processing on 1996/06/24
			// The algorithm implemented in NIH Image sets the threshold as that grey
			// value, G, for which the average of the averages of the grey values
			// below and above G is equal to G. It does this by initialising G to the
			// lowest sensible value and iterating:

			// L = the average grey value of pixels with intensities < G
			// H = the average grey value of pixels with intensities > G
			// is G = (L + H)/2?
			// yes => exit
			// no => increment G and repeat
			//
			int i, l, totl, g = 0;
			double toth, h;
			for (i = 1; i < 256; i++)
			{
				if (data[i] > 0)
				{
					g = i + 1;
					break;
				}
			}
			while (true)
			{
				l = 0;
				totl = 0;
				for (i = 0; i < g; i++)
				{
					totl = totl + data[i];
					l = l + (data[i] * i);
				}
				h = 0;
				toth = 0;
				for (i = g + 1; i < 256; i++)
				{
					toth += data[i];
					h += ((double)data[i] * i);
				}
				if (totl > 0 && toth > 0)
				{
					l /= totl;
					h /= toth;
					if (g == (int)Math.Round((l + h) / 2.0))
						break;
				}
				g++;
				if (g > 254)
					return -1;
			}
			return g;
		}

		int defaultIsoData(int[] data)
		{
			// This is the modified IsoData method used by the "Threshold" widget in "Default" mode
			int n = data.Length;
			int[] data2 = new int[n];
			int mode = 0, maxCount = 0;
			for (int i = 0; i < n; i++)
			{
				int count = data[i];
				data2[i] = data[i];
				if (data2[i] > maxCount)
				{
					maxCount = data2[i];
					mode = i;
				}
			}
			int maxCount2 = 0;
			for (int i = 0; i < n; i++)
			{
				if ((data2[i] > maxCount2) && (i != mode))
					maxCount2 = data2[i];
			}
			int hmax = maxCount;
			if ((hmax > (maxCount2 * 2)) && (maxCount2 != 0))
			{
				hmax = (int)(maxCount2 * 1.5);
				data2[mode] = hmax;
			}
			return IJIsoData(data2);
		}

		int IJIsoData(int[] data)
		{
			// This is the original ImageJ IsoData implementation, here for backward compatibility.
			int level;
			int maxValue = data.Length - 1;
			double result, sum1, sum2, sum3, sum4;
			int count0 = data[0];
			data[0] = 0; //set to zero so erased areas aren't included
			int countMax = data[maxValue];
			data[maxValue] = 0;
			int min = 0;
			while ((data[min] == 0) && (min < maxValue))
				min++;
			int max = maxValue;
			while ((data[max] == 0) && (max > 0))
				max--;
			if (min >= max)
			{
				data[0] = count0; data[maxValue] = countMax;
				level = data.Length / 2;
				return level;
			}
			int movingIndex = min;
			int inc = Math.Max(max / 40, 1);
			do
			{
				sum1 = sum2 = sum3 = sum4 = 0.0;
				for (int i = min; i <= movingIndex; i++)
				{
					sum1 += (double)i * data[i];
					sum2 += data[i];
				}
				for (int i = (movingIndex + 1); i <= max; i++)
				{
					sum3 += (double)i * data[i];
					sum4 += data[i];
				}
				result = (sum1 / sum2 + sum3 / sum4) / 2.0;
				movingIndex++;
			} while ((movingIndex + 1) <= result && movingIndex < max - 1);
			data[0] = count0; data[maxValue] = countMax;
			level = (int)Math.Round(result);
			return level;
		}


		int Li(int[] data)
		{
			// Implements Li's Minimum Cross Entropy thresholding method
			// This implementation is based on the iterative version (Ref. 2) of the algorithm.
			// 1) Li C.H. and Lee C.K. (1993) "Minimum Cross Entropy Thresholding" 
			//    Pattern Recognition, 26(4): 617-625
			// 2) Li C.H. and Tam P.K.S. (1998) "An Iterative Algorithm for Minimum 
			//    Cross Entropy Thresholding"Pattern Recognition Letters, 18(8): 771-776
			// 3) Sezgin M. and Sankur B. (2004) "Survey over Image Thresholding 
			//    Techniques and Quantitative Performance Evaluation" Journal of 
			//    Electronic Imaging, 13(1): 146-165 
			//    http://citeseer.ist.psu.edu/sezgin04survey.html
			// Ported to ImageJ plugin by G.Landini from E Celebi's fourier_0.8 routines
			int threshold;
			double num_pixels;
			double sum_back; /* sum of the background pixels at a given threshold */
			double sum_obj;  /* sum of the object pixels at a given threshold */
			double num_back; /* number of background pixels at a given threshold */
			double num_obj;  /* number of object pixels at a given threshold */
			double old_thresh;
			double new_thresh;
			double mean_back; /* mean of the background pixels at a given threshold */
			double mean_obj;  /* mean of the object pixels at a given threshold */
			double mean;  /* mean gray-level in the image */
			double tolerance; /* threshold tolerance */
			double temp;

			tolerance = 0.5;
			num_pixels = 0;
			for (int ih = 0; ih < 256; ih++)
				num_pixels += data[ih];

			/* Calculate the mean gray-level */
			mean = 0.0;
			for (int ih = 0 + 1; ih < 256; ih++) //0 + 1?
				mean += (double)ih * data[ih];
			mean /= num_pixels;
			/* Initial estimate */
			new_thresh = mean;

			do
			{
				old_thresh = new_thresh;
				threshold = (int)(old_thresh + 0.5);    /* range */
				/* Calculate the means of background and object pixels */
				/* Background */
				sum_back = 0;
				num_back = 0;
				for (int ih = 0; ih <= threshold; ih++)
				{
					sum_back += (double)ih * data[ih];
					num_back += data[ih];
				}
				mean_back = (num_back == 0 ? 0.0 : (sum_back / (double)num_back));
				/* Object */
				sum_obj = 0;
				num_obj = 0;
				for (int ih = threshold + 1; ih < 256; ih++)
				{
					sum_obj += (double)ih * data[ih];
					num_obj += data[ih];
				}
				mean_obj = (num_obj == 0 ? 0.0 : (sum_obj / (double)num_obj));

				/* Calculate the new threshold: Equation (7) in Ref. 2 */
				//new_thresh = simple_round ( ( mean_back - mean_obj ) / ( Math.log ( mean_back ) - Math.log ( mean_obj ) ) );
				//simple_round ( double x ) {
				// return ( int ) ( IS_NEG ( x ) ? x - .5 : x + .5 );
				//}
				//
				//#define IS_NEG( x ) ( ( x ) < -DBL_EPSILON ) 
				//DBL_EPSILON = 2.220446049250313E-16
				temp = (mean_back - mean_obj) / (Math.Log(mean_back) - Math.Log(mean_obj));

				if (temp < -2.220446049250313E-16)
					new_thresh = (int)(temp - 0.5);
				else
					new_thresh = (int)(temp + 0.5);
				/*  Stop the iterations when the difference between the
				new and old threshold values is less than the tolerance */
			}
			while (Math.Abs(new_thresh - old_thresh) > tolerance);
			return threshold;
		}

		int MaxEntropy(int[] data)
		{
			// Implements Kapur-Sahoo-Wong (Maximum Entropy) thresholding method
			// Kapur J.N., Sahoo P.K., and Wong A.K.C. (1985) "A New Method for
			// Gray-Level Picture Thresholding Using the Entropy of the Histogram"
			// Graphical Models and Image Processing, 29(3): 273-285
			// M. Emre Celebi
			// 06.15.2007
			// Ported to ImageJ plugin by G.Landini from E Celebi's fourier_0.8 routines
			int threshold = -1;
			int ih, it;
			int first_bin;
			int last_bin;
			double tot_ent;  /* total entropy */
			double max_ent;  /* max entropy */
			double ent_back; /* entropy of the background pixels at a given threshold */
			double ent_obj;  /* entropy of the object pixels at a given threshold */
			double[] norm_histo = new double[256]; /* normalized histogram */
			double[] P1 = new double[256]; /* cumulative normalized histogram */
			double[] P2 = new double[256];

			double total = 0;
			for (ih = 0; ih < 256; ih++)
				total += data[ih];

			for (ih = 0; ih < 256; ih++)
				norm_histo[ih] = data[ih] / total;

			P1[0] = norm_histo[0];
			P2[0] = 1.0 - P1[0];
			for (ih = 1; ih < 256; ih++)
			{
				P1[ih] = P1[ih - 1] + norm_histo[ih];
				P2[ih] = 1.0 - P1[ih];
			}

			/* Determine the first non-zero bin */
			first_bin = 0;
			for (ih = 0; ih < 256; ih++)
			{
				if (!(Math.Abs(P1[ih]) < 2.220446049250313E-16))
				{
					first_bin = ih;
					break;
				}
			}

			/* Determine the last non-zero bin */
			last_bin = 255;
			for (ih = 255; ih >= first_bin; ih--)
			{
				if (!(Math.Abs(P2[ih]) < 2.220446049250313E-16))
				{
					last_bin = ih;
					break;
				}
			}

			// Calculate the total entropy each gray-level
			// and find the threshold that maximizes it 
			max_ent = double.MinValue;

			for (it = first_bin; it <= last_bin; it++)
			{
				/* Entropy of the background pixels */
				ent_back = 0.0;
				for (ih = 0; ih <= it; ih++)
				{
					if (data[ih] != 0)
					{
						ent_back -= (norm_histo[ih] / P1[it]) * Math.Log(norm_histo[ih] / P1[it]);
					}
				}

				/* Entropy of the object pixels */
				ent_obj = 0.0;
				for (ih = it + 1; ih < 256; ih++)
				{
					if (data[ih] != 0)
					{
						ent_obj -= (norm_histo[ih] / P2[it]) * Math.Log(norm_histo[ih] / P2[it]);
					}
				}

				/* Total entropy */
				tot_ent = ent_back + ent_obj;

				// IJ.log(""+max_ent+"  "+tot_ent);
				if (max_ent < tot_ent)
				{
					max_ent = tot_ent;
					threshold = it;
				}
			}
			return threshold;
		}

		int Mean(int[] data)
		{
			// C. A. Glasbey, "An analysis of histogram-based thresholding algorithms,"
			// CVGIP: Graphical Models and Image Processing, vol. 55, pp. 532-537, 1993.
			//
			// The threshold is the mean of the greyscale data
			int threshold = -1;
			double tot = 0, sum = 0;
			for (int i = 0; i < 256; i++)
			{
				tot += data[i];
				sum += ((double)i * data[i]);
			}
			threshold = (int)Math.Floor(sum / tot);
			return threshold;
		}

		int MinErrorI(int[] data)
		{
			// Kittler and J. Illingworth, "Minimum error thresholding," Pattern Recognition, vol. 19, pp. 41-47, 1986.
			// C. A. Glasbey, "An analysis of histogram-based thresholding algorithms," CVGIP: Graphical Models and Image Processing, vol. 55, pp. 532-537, 1993.
			// Ported to ImageJ plugin by G.Landini from Antti Niemisto's Matlab code (GPL)
			// Original Matlab code Copyright (C) 2004 Antti Niemisto
			// See http://www.cs.tut.fi/~ant/histthresh/ for an excellent slide presentation
			// and the original Matlab code.

			int threshold = Mean(data); //Initial estimate for the threshold is found with the MEAN algorithm.
			int Tprev = -2;
			double mu, nu, p, q, sigma2, tau2, w0, w1, w2, sqterm, temp;
			//int counter=1;
			while (threshold != Tprev)
			{
				//Calculate some statistics.
				mu = B(data, threshold) / A(data, threshold);
				nu = (B(data, data.Length - 1) - B(data, threshold)) / (A(data, data.Length - 1) - A(data, threshold));
				p = A(data, threshold) / A(data, data.Length - 1);
				q = (A(data, data.Length - 1) - A(data, threshold)) / A(data, data.Length - 1);
				sigma2 = C(data, threshold) / A(data, threshold) - (mu * mu);
				tau2 = (C(data, data.Length - 1) - C(data, threshold)) / (A(data, data.Length - 1) - A(data, threshold)) - (nu * nu);

				//The terms of the quadratic equation to be solved.
				w0 = 1.0 / sigma2 - 1.0 / tau2;
				w1 = mu / sigma2 - nu / tau2;
				w2 = (mu * mu) / sigma2 - (nu * nu) / tau2 + Math.Log10((sigma2 * (q * q)) / (tau2 * (p * p)));

				//If the next threshold would be imaginary, return with the current one.
				sqterm = (w1 * w1) - w0 * w2;
				if (sqterm < 0)
				{
					IJ.log("MinError(I): not converging.");
					return threshold;
				}

				//The updated threshold is the integer part of the solution of the quadratic equation.
				Tprev = threshold;
				temp = (w1 + Math.Sqrt(sqterm)) / w0;

				if (double.IsNaN(temp))
					threshold = Tprev;
				else
					threshold = (int)Math.Floor(temp);
			}
			return threshold;
		}

		private double A(int[] y, int j)
		{
			if (j >= y.Length) j = y.Length - 1;
			double x = 0;
			for (int i = 0; i <= j; i++)
				x += y[i];
			return x;
		}

		private double B(int[] y, int j)
		{
			if (j >= y.Length) j = y.Length - 1;
			double x = 0;
			for (int i = 0; i <= j; i++)
				x += i * y[i];
			return x;
		}

		private double C(int[] y, int j)
		{
			if (j >= y.Length) j = y.Length - 1;
			double x = 0;
			for (int i = 0; i <= j; i++)
				x += i * i * y[i];
			return x;
		}

		int Minimum(int[] data)
		{
			// J. M. S. Prewitt and M. L. Mendelsohn, "The analysis of cell images," in
			// Annals of the New York Academy of Sciences, vol. 128, pp. 1035-1053, 1966.
			// ported to ImageJ plugin by G.Landini from Antti Niemisto's Matlab code (GPL)
			// Original Matlab code Copyright (C) 2004 Antti Niemisto
			// See http://www.cs.tut.fi/~ant/histthresh/ for an excellent slide presentation
			// and the original Matlab code.
			//
			// Assumes a bimodal histogram. The histogram needs is smoothed (using a
			// running average of size 3, iteratively) until there are only two local maxima.
			// Threshold t is such that yt-1 > yt <= yt+1.
			// Images with histograms having extremely unequal peaks or a broad and
			// flat valleys are unsuitable for this method.
			int iter = 0;
			int threshold = -1;
			double[] iHisto = new double[256];
			for (int i = 0; i < 256; i++)
				iHisto[i] = (double)data[i];
			double[] tHisto = new double[iHisto.Length];

			while (!bimodalTest(iHisto))
			{
				//smooth with a 3 point running mean filter
				for (int i = 1; i < 255; i++)
					tHisto[i] = (iHisto[i - 1] + iHisto[i] + iHisto[i + 1]) / 3;
				tHisto[0] = (iHisto[0] + iHisto[1]) / 3; //0 outside
				tHisto[255] = (iHisto[254] + iHisto[255]) / 3; //0 outside
				Array.Copy(tHisto, 0, iHisto, 0, iHisto.Length);
				iter++;
				if (iter > 10000)
				{
					threshold = -1;
					IJ.log("Minimum: threshold not found after 10000 iterations.");
					return threshold;
				}
			}
			// The threshold is the minimum between the two peaks.
			for (int i = 1; i < 255; i++)
			{
				if (iHisto[i - 1] > iHisto[i] && iHisto[i + 1] >= iHisto[i])
				{
					threshold = i;
					break;
				}
			}
			return threshold;
		}

		int Moments(int[] data)
		{
			//  W. Tsai, "Moment-preserving thresholding: a new approach," Computer Vision,
			// Graphics, and Image Processing, vol. 29, pp. 377-393, 1985.
			// Ported to ImageJ plugin by G.Landini from the the open source project FOURIER 0.8
			// by  M. Emre Celebi , Department of Computer Science,  Louisiana State University in Shreveport
			// Shreveport, LA 71115, USA
			//  http://sourceforge.net/projects/fourier-ipal
			//  http://www.lsus.edu/faculty/~ecelebi/fourier.htm
			double total = 0;
			double m0 = 1.0, m1 = 0.0, m2 = 0.0, m3 = 0.0, sum = 0.0, p0 = 0.0;
			double cd, c0, c1, z0, z1;  /* auxiliary variables */
			int threshold = -1;

			double[] histo = new double[256];

			for (int i = 0; i < 256; i++)
				total += data[i];

			for (int i = 0; i < 256; i++)
				histo[i] = (double)(data[i] / total); //normalised histogram

			/* Calculate the first, second, and third order moments */
			for (int i = 0; i < 256; i++)
			{
				double di = i;
				m1 += di * histo[i];
				m2 += di * di * histo[i];
				m3 += di * di * di * histo[i];
			}
			/* 
			First 4 moments of the gray-level image should match the first 4 moments
			of the target binary image. This leads to 4 equalities whose solutions 
			are given in the Appendix of Ref. 1 
			*/
			cd = m0 * m2 - m1 * m1;
			c0 = (-m2 * m2 + m1 * m3) / cd;
			c1 = (m0 * -m3 + m2 * m1) / cd;
			z0 = 0.5 * (-c1 - Math.Sqrt(c1 * c1 - 4.0 * c0));
			z1 = 0.5 * (-c1 + Math.Sqrt(c1 * c1 - 4.0 * c0));
			p0 = (z1 - m1) / (z1 - z0);  /* Fraction of the object pixels in the target binary image */

			// The threshold is the gray-level closest  
			// to the p0-tile of the normalized histogram 
			sum = 0;
			for (int i = 0; i < 256; i++)
			{
				sum += histo[i];
				if (sum > p0)
				{
					threshold = i;
					break;
				}
			}
			return threshold;
		}

		int Otsu(int[] data)
		{
			// Otsu's threshold algorithm
			// C++ code by Jordan Bevik <Jordan.Bevic@qtiworld.com>
			// ported to ImageJ plugin by G.Landini
			int k, kStar;  // k = the current threshold; kStar = optimal threshold
			double N1, N;    // N1 = # points with intensity <=k; N = total number of points
			double BCV, BCVmax; // The current Between Class Variance and maximum BCV
			double num, denom;  // temporary bookeeping
			double Sk;  // The total intensity for all histogram points <=k
			double S, L = 256; // The total intensity of the image

			// Initialize values:
			S = N = 0;
			for (k = 0; k < L; k++)
			{
				S += (double)k * data[k];   // Total histogram intensity
				N += data[k];       // Total number of data points
			}

			Sk = 0;
			N1 = data[0]; // The entry for zero intensity
			BCV = 0;
			BCVmax = 0;
			kStar = 0;

			// Look at each possible threshold value,
			// calculate the between-class variance, and decide if it's a max
			for (k = 1; k < L - 1; k++)
			{ // No need to check endpoints k = 0 or k = L-1
				Sk += (double)k * data[k];
				N1 += data[k];

				// The float casting here is to avoid compiler warning about loss of precision and
				// will prevent overflow in the case of large saturated images
				denom = (double)(N1) * (N - N1); // Maximum value of denom is (N^2)/4 =  approx. 3E10

				if (denom != 0)
				{
					// Float here is to avoid loss of precision when dividing
					num = ((double)N1 / N) * S - Sk;    // Maximum value of num =  255*N = approx 8E7
					BCV = (num * num) / denom;
				}
				else
					BCV = 0;

				if (BCV >= BCVmax)
				{ // Assign the best threshold found so far
					BCVmax = BCV;
					kStar = k;
				}
			}
			// kStar += 1;	// Use QTI convention that intensity -> 1 if intensity >= k
			// (the algorithm was developed for I-> 1 if I <= k.)
			return kStar;
		}


		int Percentile(int[] data)
		{
			// W. Doyle, "Operation useful for similarity-invariant pattern recognition,"
			// Journal of the Association for Computing Machinery, vol. 9,pp. 259-267, 1962.
			// ported to ImageJ plugin by G.Landini from Antti Niemisto's Matlab code (GPL)
			// Original Matlab code Copyright (C) 2004 Antti Niemisto
			// See http://www.cs.tut.fi/~ant/histthresh/ for an excellent slide presentation
			// and the original Matlab code.

			int iter = 0;
			int threshold = -1;
			double ptile = 0.5; // default fraction of foreground pixels
			double[] avec = new double[256];

			for (int i = 0; i < 256; i++)
				avec[i] = 0.0;

			double total = partialSum(data, 255);
			double temp = 1.0;
			for (int i = 0; i < 256; i++)
			{
				avec[i] = Math.Abs((partialSum(data, i) / total) - ptile);
				//IJ.log("Ptile["+i+"]:"+ avec[i]);
				if (avec[i] < temp)
				{
					temp = avec[i];
					threshold = i;
				}
			}
			return threshold;
		}


		double partialSum(int[] y, int j)
		{
			double x = 0;
			for (int i = 0; i <= j; i++)
				x += y[i];
			return x;
		}


		int RenyiEntropy(int[] data)
		{
			// Kapur J.N., Sahoo P.K., and Wong A.K.C. (1985) "A New Method for
			// Gray-Level Picture Thresholding Using the Entropy of the Histogram"
			// Graphical Models and Image Processing, 29(3): 273-285
			// M. Emre Celebi
			// 06.15.2007
			// Ported to ImageJ plugin by G.Landini from E Celebi's fourier_0.8 routines

			int threshold;
			int opt_threshold;

			int ih, it;
			int first_bin;
			int last_bin;
			int tmp_var;
			int t_star1, t_star2, t_star3;
			int beta1, beta2, beta3;
			double alpha;/* alpha parameter of the method */
			double term;
			double tot_ent;  /* total entropy */
			double max_ent;  /* max entropy */
			double ent_back; /* entropy of the background pixels at a given threshold */
			double ent_obj;  /* entropy of the object pixels at a given threshold */
			double omega;
			double[] norm_histo = new double[256]; /* normalized histogram */
			double[] P1 = new double[256]; /* cumulative normalized histogram */
			double[] P2 = new double[256];

			double total = 0;
			for (ih = 0; ih < 256; ih++)
				total += data[ih];

			for (ih = 0; ih < 256; ih++)
				norm_histo[ih] = data[ih] / total;

			P1[0] = norm_histo[0];
			P2[0] = 1.0 - P1[0];
			for (ih = 1; ih < 256; ih++)
			{
				P1[ih] = P1[ih - 1] + norm_histo[ih];
				P2[ih] = 1.0 - P1[ih];
			}

			/* Determine the first non-zero bin */
			first_bin = 0;
			for (ih = 0; ih < 256; ih++)
			{
				if (!(Math.Abs(P1[ih]) < 2.220446049250313E-16))
				{
					first_bin = ih;
					break;
				}
			}

			/* Determine the last non-zero bin */
			last_bin = 255;
			for (ih = 255; ih >= first_bin; ih--)
			{
				if (!(Math.Abs(P2[ih]) < 2.220446049250313E-16))
				{
					last_bin = ih;
					break;
				}
			}

			/* Maximum Entropy Thresholding - BEGIN */
			/* ALPHA = 1.0 */
			/* Calculate the total entropy each gray-level
			and find the threshold that maximizes it 
			*/
			threshold = 0; // was MIN_INT in original code, but if an empty image is processed it gives an error later on.
			max_ent = 0.0;

			for (it = first_bin; it <= last_bin; it++)
			{
				/* Entropy of the background pixels */
				ent_back = 0.0;
				for (ih = 0; ih <= it; ih++)
				{
					if (data[ih] != 0)
					{
						ent_back -= (norm_histo[ih] / P1[it]) * Math.Log(norm_histo[ih] / P1[it]);
					}
				}

				/* Entropy of the object pixels */
				ent_obj = 0.0;
				for (ih = it + 1; ih < 256; ih++)
				{
					if (data[ih] != 0)
					{
						ent_obj -= (norm_histo[ih] / P2[it]) * Math.Log(norm_histo[ih] / P2[it]);
					}
				}

				/* Total entropy */
				tot_ent = ent_back + ent_obj;

				// IJ.log(""+max_ent+"  "+tot_ent);

				if (max_ent < tot_ent)
				{
					max_ent = tot_ent;
					threshold = it;
				}
			}
			t_star2 = threshold;

			/* Maximum Entropy Thresholding - END */
			threshold = 0; //was MIN_INT in original code, but if an empty image is processed it gives an error later on.
			max_ent = 0.0;
			alpha = 0.5;
			term = 1.0 / (1.0 - alpha);
			for (it = first_bin; it <= last_bin; it++)
			{
				/* Entropy of the background pixels */
				ent_back = 0.0;
				for (ih = 0; ih <= it; ih++)
					ent_back += Math.Sqrt(norm_histo[ih] / P1[it]);

				/* Entropy of the object pixels */
				ent_obj = 0.0;
				for (ih = it + 1; ih < 256; ih++)
					ent_obj += Math.Sqrt(norm_histo[ih] / P2[it]);

				/* Total entropy */
				tot_ent = term * ((ent_back * ent_obj) > 0.0 ? Math.Log(ent_back * ent_obj) : 0.0);

				if (tot_ent > max_ent)
				{
					max_ent = tot_ent;
					threshold = it;
				}
			}

			t_star1 = threshold;

			threshold = 0; //was MIN_INT in original code, but if an empty image is processed it gives an error later on.
			max_ent = 0.0;
			alpha = 2.0;
			term = 1.0 / (1.0 - alpha);
			for (it = first_bin; it <= last_bin; it++)
			{
				/* Entropy of the background pixels */
				ent_back = 0.0;
				for (ih = 0; ih <= it; ih++)
					ent_back += (norm_histo[ih] * norm_histo[ih]) / (P1[it] * P1[it]);

				/* Entropy of the object pixels */
				ent_obj = 0.0;
				for (ih = it + 1; ih < 256; ih++)
					ent_obj += (norm_histo[ih] * norm_histo[ih]) / (P2[it] * P2[it]);

				/* Total entropy */
				tot_ent = term * ((ent_back * ent_obj) > 0.0 ? Math.Log(ent_back * ent_obj) : 0.0);

				if (tot_ent > max_ent)
				{
					max_ent = tot_ent;
					threshold = it;
				}
			}

			t_star3 = threshold;

			/* Sort t_star values */
			if (t_star2 < t_star1)
			{
				tmp_var = t_star1;
				t_star1 = t_star2;
				t_star2 = tmp_var;
			}
			if (t_star3 < t_star2)
			{
				tmp_var = t_star2;
				t_star2 = t_star3;
				t_star3 = tmp_var;
			}
			if (t_star2 < t_star1)
			{
				tmp_var = t_star1;
				t_star1 = t_star2;
				t_star2 = tmp_var;
			}

			/* Adjust beta values */
			if (Math.Abs(t_star1 - t_star2) <= 5)
			{
				if (Math.Abs(t_star2 - t_star3) <= 5)
				{
					beta1 = 1;
					beta2 = 2;
					beta3 = 1;
				}
				else
				{
					beta1 = 0;
					beta2 = 1;
					beta3 = 3;
				}
			}
			else
			{
				if (Math.Abs(t_star2 - t_star3) <= 5)
				{
					beta1 = 3;
					beta2 = 1;
					beta3 = 0;
				}
				else
				{
					beta1 = 1;
					beta2 = 2;
					beta3 = 1;
				}
			}
			//IJ.log(""+t_star1+" "+t_star2+" "+t_star3);
			/* Determine the optimal threshold value */
			omega = P1[t_star3] - P1[t_star1];
			opt_threshold = (int)(t_star1 * (P1[t_star1] + 0.25 * omega * beta1) + 0.25 * t_star2 * omega * beta2 + t_star3 * (P2[t_star3] + 0.25 * omega * beta3));

			return opt_threshold;
		}


		int Shanbhag(int[] data)
		{
			// Shanhbag A.G. (1994) "Utilization of Information Measure as a Means of
			//  Image Thresholding" Graphical Models and Image Processing, 56(5): 414-419
			// Ported to ImageJ plugin by G.Landini from E Celebi's fourier_0.8 routines
			int threshold;
			int ih, it;
			int first_bin;
			int last_bin;
			double term;
			double tot_ent;  /* total entropy */
			double min_ent;  /* max entropy */
			double ent_back; /* entropy of the background pixels at a given threshold */
			double ent_obj;  /* entropy of the object pixels at a given threshold */
			double[] norm_histo = new double[256]; /* normalized histogram */
			double[] P1 = new double[256]; /* cumulative normalized histogram */
			double[] P2 = new double[256];

			double total = 0;
			for (ih = 0; ih < 256; ih++)
				total += data[ih];

			for (ih = 0; ih < 256; ih++)
				norm_histo[ih] = data[ih] / total;

			P1[0] = norm_histo[0];
			P2[0] = 1.0 - P1[0];
			for (ih = 1; ih < 256; ih++)
			{
				P1[ih] = P1[ih - 1] + norm_histo[ih];
				P2[ih] = 1.0 - P1[ih];
			}

			/* Determine the first non-zero bin */
			first_bin = 0;
			for (ih = 0; ih < 256; ih++)
			{
				if (!(Math.Abs(P1[ih]) < 2.220446049250313E-16))
				{
					first_bin = ih;
					break;
				}
			}

			/* Determine the last non-zero bin */
			last_bin = 255;
			for (ih = 255; ih >= first_bin; ih--)
			{
				if (!(Math.Abs(P2[ih]) < 2.220446049250313E-16))
				{
					last_bin = ih;
					break;
				}
			}

			// Calculate the total entropy each gray-level
			// and find the threshold that maximizes it 
			threshold = -1;
			min_ent = double.MaxValue;

			for (it = first_bin; it <= last_bin; it++)
			{
				/* Entropy of the background pixels */
				ent_back = 0.0;
				term = 0.5 / P1[it];
				for (ih = 1; ih <= it; ih++)
				{ //0+1?
					ent_back -= norm_histo[ih] * Math.Log(1.0 - term * P1[ih - 1]);
				}
				ent_back *= term;

				/* Entropy of the object pixels */
				ent_obj = 0.0;
				term = 0.5 / P2[it];
				for (ih = it + 1; ih < 256; ih++)
				{
					ent_obj -= norm_histo[ih] * Math.Log(1.0 - term * P2[ih]);
				}
				ent_obj *= term;

				/* Total entropy */
				tot_ent = Math.Abs(ent_back - ent_obj);

				if (tot_ent < min_ent)
				{
					min_ent = tot_ent;
					threshold = it;
				}
			}
			return threshold;
		}


		int Triangle(int[] data)
		{
			//  Zack, G. W., Rogers, W. E. and Latt, S. A., 1977,
			//  Automatic Measurement of Sister Chromatid Exchange Frequency,
			// Journal of Histochemistry and Cytochemistry 25 (7), pp. 741-753
			//
			//  modified from Johannes Schindelin plugin
			// 
			// find min and max
			int min = 0, dmax = 0, max = 0, min2 = 0;
			for (int i = 0; i < data.Length; i++)
			{
				if (data[i] > 0)
				{
					min = i;
					break;
				}
			}
			if (min > 0) min--; // line to the (p==0) point, not to data[min]

			// The Triangle algorithm cannot tell whether the data is skewed to one side or another.
			// This causes a problem as there are 2 possible thresholds between the max and the 2 extremes
			// of the histogram.
			// Here I propose to find out to which side of the max point the data is furthest, and use that as
			//  the other extreme.
			for (int i = 255; i > 0; i--)
			{
				if (data[i] > 0)
				{
					min2 = i;
					break;
				}
			}
			if (min2 < 255) min2++; // line to the (p==0) point, not to data[min]

			for (int i = 0; i < 256; i++)
			{
				if (data[i] > dmax)
				{
					max = i;
					dmax = data[i];
				}
			}
			// find which is the furthest side
			//IJ.log(""+min+" "+max+" "+min2);
			bool inverted = false;
			if ((max - min) < (min2 - max))
			{
				// reverse the histogram
				//IJ.log("Reversing histogram.");
				inverted = true;
				int left = 0;          // index of leftmost element
				int right = 255; // index of rightmost element
				while (left < right)
				{
					// exchange the left and right elements
					int temp = data[left];
					data[left] = data[right];
					data[right] = temp;
					// move the bounds toward the center
					left++;
					right--;
				}
				min = 255 - min2;
				max = 255 - max;
			}

			if (min == max)
			{
				//IJ.log("Triangle:  min == max.");
				return min;
			}

			// describe line by nx * x + ny * y - d = 0
			double nx, ny, d;
			// nx is just the max frequency as the other point has freq=0
			nx = data[max];   //-min; // data[min]; //  lowest value bmin = (p=0)% in the image
			ny = min - max;
			d = Math.Sqrt(nx * nx + ny * ny);
			nx /= d;
			ny /= d;
			d = nx * min + ny * data[min];

			// find split point
			int split = min;
			double splitDistance = 0;
			for (int i = min + 1; i <= max; i++)
			{
				double newDistance = nx * i + ny * data[i] - d;
				if (newDistance > splitDistance)
				{
					split = i;
					splitDistance = newDistance;
				}
			}
			split--;

			if (inverted)
			{
				// The histogram might be used for something else, so let's reverse it back
				int left = 0;
				int right = 255;
				while (left < right)
				{
					int temp = data[left];
					data[left] = data[right];
					data[right] = temp;
					left++;
					right--;
				}
				return (255 - split);
			}
			else
				return split;
		}


		int Yen(int[] data)
		{
			// Implements Yen  thresholding method
			// 1) Yen J.C., Chang F.J., and Chang S. (1995) "A New Criterion 
			//    for Automatic Multilevel Thresholding" IEEE Trans. on Image 
			//    Processing, 4(3): 370-378
			// 2) Sezgin M. and Sankur B. (2004) "Survey over Image Thresholding 
			//    Techniques and Quantitative Performance Evaluation" Journal of 
			//    Electronic Imaging, 13(1): 146-165
			//    http://citeseer.ist.psu.edu/sezgin04survey.html
			//
			// M. Emre Celebi
			// 06.15.2007
			// Ported to ImageJ plugin by G.Landini from E Celebi's fourier_0.8 routines
			int threshold;
			int ih, it;
			double crit;
			double max_crit;
			double[] norm_histo = new double[256]; /* normalized histogram */
			double[] P1 = new double[256]; /* cumulative normalized histogram */
			double[] P1_sq = new double[256];
			double[] P2_sq = new double[256];

			double total = 0;
			for (ih = 0; ih < 256; ih++)
				total += data[ih];

			for (ih = 0; ih < 256; ih++)
				norm_histo[ih] = data[ih] / total;

			P1[0] = norm_histo[0];
			for (ih = 1; ih < 256; ih++)
				P1[ih] = P1[ih - 1] + norm_histo[ih];

			P1_sq[0] = norm_histo[0] * norm_histo[0];
			for (ih = 1; ih < 256; ih++)
				P1_sq[ih] = P1_sq[ih - 1] + norm_histo[ih] * norm_histo[ih];

			P2_sq[255] = 0.0;
			for (ih = 254; ih >= 0; ih--)
				P2_sq[ih] = P2_sq[ih + 1] + norm_histo[ih + 1] * norm_histo[ih + 1];

			/* Find the threshold that maximizes the criterion */
			threshold = -1;
			max_crit = double.MinValue;
			for (it = 0; it < 256; it++)
			{
				crit = -1.0 * ((P1_sq[it] * P2_sq[it]) > 0.0 ? Math.Log(P1_sq[it] * P2_sq[it]) : 0.0) + 2 * ((P1[it] * (1.0 - P1[it])) > 0.0 ? Math.Log(P1[it] * (1.0 - P1[it])) : 0.0);
				if (crit > max_crit)
				{
					max_crit = crit;
					threshold = it;
				}
			}
			return threshold;
		}
	}
}
